using System;
using System.Runtime.Serialization;
using Aima.Core.util;


namespace Aima.Core.Util.math
{
/**
 * LU Decomposition.
 * For an m-by-n matrix A with m >= n, the LU decomposition is an m-by-n unit
 * lower triangular matrix L, an n-by-n upper triangular matrix U, and a
 * permutation vector piv of length m so that A(piv,:) = L*U. 
 * If m < n, then L  * is m-by-m and U is m-by-n.
 * used is in the solution of square systems of simultaneous linear
 * equations. This will fail if isNonsingular() returns false.
 */

    public class LuDecomposition : ISerializable
    {
        // internal array storage.
        private readonly double[][] _lu;

        // Row and column dimensions, pivot vector, and pivot sign.
        private readonly int _m;
        private readonly int _n;
        private readonly int[] _piv;
        private readonly int _pivsign;

        public  int M {get { return _m; }}
        public  int N {get { return _n; }}
        public int[] Piv {get { return _piv; }}
        public int PivSign {get { return _pivsign; }}
     
        public LuDecomposition(Matrix matrix)
        {
            // Use a "left-looking", dot-product, Crout/Doolittle algorithm.
            _lu = matrix.GetArrayCopy();
            _m = matrix.RowDimension;
            _n = matrix.ColumnDimension;
            _piv = new int[_m];
            for (int i = 0; i < _m; i++)
            {
                _piv[i] = i;
            }
            _pivsign = 1;
            var lUcolj = new double[_m];

            for (int j = 0; j < _n; j++)
            {
                // Make a copy of the j-th column to localize references.
                for (int i = 0; i < _m; i++)
                {
                    lUcolj[i] = _lu[i][j];
                }

                // Apply previous transformations.

                for (int i = 0; i < _m; i++)
                {
                    double[] lUrowi = _lu[i];

                    // Most of the time is spent in the following dot product.

                    int kmax = Math.Min(i, j);
                    double s = 0.0;
                    for (int k = 0; k < kmax; k++)
                    {
                        s += lUrowi[k]*lUcolj[k];
                    }

                    lUrowi[j] = lUcolj[i] -= s;
                }

                // Find pivot and exchange if necessary.

                int p = j;
                for (int i = j + 1; i < _m; i++)
                {
                    if (Math.Abs(lUcolj[i]) > Math.Abs(lUcolj[p]))
                    {
                        p = i;
                    }
                }
                if (p != j)
                {
                    for (int k = 0; k < _n; k++)
                    {
                        double t = _lu[p][k];
                        _lu[p][k] = _lu[j][k];
                        _lu[j][k] = t;
                    }
                    int k2 = _piv[p];
                    _piv[p] = _piv[j];
                    _piv[j] = k2;
                    _pivsign = -_pivsign;
                }

                // Compute multipliers.
                // TODO: fix float comparison to use an epsilon
                if (j < _m & Math.Abs(_lu[j][j] - 0.0) > Consts.Epsilon)
                {
                    for (int i = j + 1; i < _m; i++)
                    {
                        _lu[i][j] /= _lu[j][j];
                    }
                }
            }
        }

        #region ISerializable 

        public void GetObjectData(SerializationInfo info, StreamingContext context)
        {
            info.AddValue("M", M);
            info.AddValue("N", N);
            info.AddValue("Piv", Piv);
            info.AddValue("PivSign", PivSign);
        }

        #endregion

        /*
	 * ------------------------ Temporary, experimental code.
	 * ------------------------\
	 * 
	 * \ LU Decomposition, computed by Gaussian elimination. <P> This
	 * constructor computes L and U with the "daxpy"-based elimination algorithm
	 * used in LINPACK and MATLAB. In Java, we suspect the dot-product, Crout
	 * algorithm will be faster. We have temporarily included this constructor
	 * until timing experiments confirm this suspicion. <P> @param A Rectangular
	 * matrix @param linpackflag Use Gaussian elimination. Actual value ignored.
	 * 
	 * @return Structure to access L, U and piv. \
	 * 
	 * public LUDecomposition (Matrix A, int linpackflag) { // Initialize. LU =
	 * A.getArrayCopy(); m = A.getRowDimension(); n = A.getColumnDimension();
	 * piv = new int[m]; for (int i = 0; i < m; i++) { piv[i] = i; } pivsign =
	 * 1; // Main loop. for (int k = 0; k < n; k++) { // Find pivot. int p = k;
	 * for (int i = k+1; i < m; i++) { if (Math.abs(LU[i][k]) >
	 * Math.abs(LU[p][k])) { p = i; } } // Exchange if necessary. if (p != k) {
	 * for (int j = 0; j < n; j++) { double t = LU[p][j]; LU[p][j] = LU[k][j];
	 * LU[k][j] = t; } int t = piv[p]; piv[p] = piv[k]; piv[k] = t; pivsign =
	 * -pivsign; } // Compute multipliers and eliminate k-th column. if
	 * (LU[k][k] != 0.0) { for (int i = k+1; i < m; i++) { LU[i][k] /= LU[k][k];
	 * for (int j = k+1; j < n; j++) { LU[i][j] -= LU[i][k]LU[k][j]; } } } } } \
	 * ------------------------ End of temporary code. ------------------------
	 */

        /*
	 * ------------------------ Public Methods ------------------------
	 */

        // Is the matrix nonsingular? return true if U, and hence A, is nonsingular.
	    public bool IsNonsingular()
	    {
	        
            for (int j = 0; j < _n; j++)
            {
                if (Math.Abs(_lu[j][j] - 0) < Consts.Epsilon)
                    return false;
            }
            return true;
        }

        // Return lower triangular factor L
        public Matrix GetL()
        {
            var x = new Matrix(_m, _n);
            double[][] l = x.MatrixArray;
            for (int i = 0; i < _m; i++)
            {
                for (int j = 0; j < _n; j++)
                {
                    if (i > j)
                    {
                        l[i][j] = _lu[i][j];
                    }
                    else if (i == j)
                    {
                        l[i][j] = 1.0;
                    }
                    else
                    {
                        l[i][j] = 0.0;
                    }
                }
            }
            return x;
        }

        // Return upper triangular factor U
        public Matrix GetU()
        {
            var x = new Matrix(_n, _n);
            double[][] u = x.MatrixArray;
            for (int i = 0; i < _n; i++)
            {
                for (int j = 0; j < _n; j++)
                {
                    if (i <= j)
                    {
                        u[i][j] = _lu[i][j];
                    }
                    else
                    {
                        u[i][j] = 0.0;
                    }
                }
            }
            return x;
        }

        // Return pivot permutation vector piv
	    public int[] GetPivot()
        {
            var p = new int[_m];
            for (int i = 0; i < _m; i++)
            {
                p[i] = _piv[i];
            }
            return p;
        }

        // Return pivot permutation vector as a one-dimensional double array
	    public double[] GetDoublePivot()
        {
            var vals = new double[_m];
            for (int i = 0; i < _m; i++)
            {
                vals[i] = _piv[i];
            }
            return vals;
        }

        // Determinant
	    public double Det()
        {
            if (_m != _n)
            {
                throw new ArgumentException("Matrix must be square.");
            }
            double d = _pivsign;
            for (int j = 0; j < _n; j++)
            {
                d *= _lu[j][j];
            }
            return d;
        }

        // Solve A*X = B & return X so that L*U*X = B(piv,:)
        public Matrix Solve(Matrix b)
        {
            if (b.RowDimension != _m)
            {
                throw new ArgumentException("Matrix row dimensions must agree.");
            }
            if (!IsNonsingular())
            {
                throw new Exception("Matrix is singular.");
            }

            // Copy right hand side with pivoting
            int nx = b.ColumnDimension;
            Matrix xmat = b.GetMatrix(_piv, 0, nx - 1);
            double[][] x = xmat.MatrixArray;

            // Solve L*Y = B(piv,:)
            for (int k = 0; k < _n; k++)
            {
                for (int i = k + 1; i < _n; i++)
                {
                    for (int j = 0; j < nx; j++)
                    {
                        x[i][j] -= x[k][j]*_lu[i][k];
                    }
                }
            }
            // Solve U*X = Y;
            for (int k = _n - 1; k >= 0; k--)
            {
                for (int j = 0; j < nx; j++)
                {
                    x[k][j] /= _lu[k][k];
                }
                for (int i = 0; i < k; i++)
                {
                    for (int j = 0; j < nx; j++)
                    {
                        x[i][j] -= x[k][j]*_lu[i][k];
                    }
                }
            }
            return xmat;
        }
    }
}